#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#include <string.h>
#include <unistd.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#define N 5
#define M 50

void histograma(gsl_vector *eval, gsl_matrix *hist, int n);
void obterDiscos(gsl_matrix* A, gsl_matrix* info);
void imprimirDiscos(gsl_matrix* info, int n);
void histDiscos(gsl_matrix* info, gsl_matrix* histD, int m);

int main(void)
{
	
	int i,j;
	gsl_matrix *A, *info;
	gsl_matrix *hist;
	double tmp;
	FILE *fp;
	char str[100];
	int m = 100;

	A = gsl_matrix_alloc(N, N);
	hist = gsl_matrix_alloc(m,2);
	info = gsl_matrix_alloc(N,2);

	/* Define a matriz A */
	for ( i = 0 ; i < N ; i++ )
	{
		for ( j = 0 ; j < N ; j++ )
		{
			gsl_matrix_set(A, i, j, 1.0/(i + j + 1) );
		}
	}
	
	/* Calcular os autovalores desta matriz */
	
	// gsl_eigen_symm(A, eval, w);
	
	// gsl_sort_vector_smallest (dest, N/frac, eval);
	
/*	for ( i = 0 ; i < N/frac ; i++ )
	{
		gsl_vector_set(subset, i, dest[i]);
	}*/
	
	/* Fazer um histograma dos autovalores */
//	histograma(subset, hist, N/frac);

	/* Imprimir o resultado */
// 	printf("Autovalores\n");
// 	for ( i = 0 ; i < N ; i++ )
// 	{
// 		printf("%18.8le\n", gsl_vector_get(eval, i));
// 	}

	//printf("Histograma\n");
	
/*	sprintf(str, "Histograma_Autovalores");
	fp = fopen(str, "w+");
	
	for ( i = 0 ; i < M ; i++ )
	{
		fprintf(fp, "%-18.8le%4d\n", gsl_matrix_get(hist, i, 0), (int)gsl_matrix_get(hist, i, 1));
	}
	fclose(fp);

	gsl_eigen_symm_free(w);*/
	
	/* Tomar a matriz e calcular o centro e os raios dos discos */
	
	obterDiscos(A, info);
//	imprimirDiscos(info, N);
	/* Organizar a informação num histograma */
	histDiscos(info, hist, 100);
	imprimirDiscos(hist, 100);

	
	return 0;

}

void histDiscos(gsl_matrix* info, gsl_matrix* hist, int m)
{
	double max, min;
	double tmp;
	double interval, points;
	double val[m];
	int i, j;
	
	
	max = gsl_matrix_get(info, 0, 0) + gsl_matrix_get(info, 0, 1);
	min = gsl_matrix_get(info, 0, 0) - gsl_matrix_get(info, 0, 1);
	
	/* Valor maximo e minimo de todo o intervalo */
	for ( i = 1 ; i < N ; i++ )
	{
		tmp = gsl_matrix_get(info, i, 0) + gsl_matrix_get(info, i, 1); 
		if ( tmp > max)
		{
			max = tmp; 
		}
		tmp = gsl_matrix_get(info, i, 0) - gsl_matrix_get(info, i, 1); 
		if ( tmp < min)
		{
			min = tmp; 
		}
	}
	
	interval = fabs(max - min);
	points = interval/(1.0*m);
	
	/* Vetor val[i] contém os sub-intervalos */
	for ( i = 0 ; i < m ; i++ )
	{
		val[i] = min + i*points;
	}
	
	gsl_matrix_set_zero(hist);
	
	/* Para cada centro e raio, gerar valores máximos e mínimo */
	/* Somar 1 no histograma se min < val[i] < max */
	for ( i = 0 ; i < N ; i++ )
	{
		max = gsl_matrix_get(info, i, 0) + gsl_matrix_get(info, i, 1);
		min = gsl_matrix_get(info, i, 0) - gsl_matrix_get(info, i, 1);
		for ( j = 0 ; j < m ; j++ )
		{
			if ( (val[j] >= min) && (val[j] <= max) )
			{
				gsl_matrix_set(hist, j, 0, val[j]);
				gsl_matrix_set(hist, j, 1, gsl_matrix_get(hist, j, 1) + 1);
			}
		}
	}
		
}

void obterDiscos(gsl_matrix* A, gsl_matrix* info)
{
	int i,j;
	double cont;
	for (i = 0 ; i < N ; i++ )
	{
		gsl_matrix_set(info, i, 0, gsl_matrix_get(A, i, i));
		cont = 0;
		for ( j = 0 ; j < N ; j++ )
		{
			if (j == i) continue;
			cont = cont + fabs(gsl_matrix_get(A, i, j));
		}
		gsl_matrix_set(info, i, 1, cont);
	}
}

void imprimirDiscos(gsl_matrix* info, int n)
{
	int i;
	for ( i = 0 ; i < n ; i++ )
	{
		printf("%18.8le%18.8le\n", gsl_matrix_get(info, i, 0), gsl_matrix_get(info, i, 1));
	}
}

// void histograma(gsl_vector *eval, gsl_matrix *hist, int n)
// {
// 	double max, min;
// 	double interval, points, val1, val2, valm;
// 	int i, j, cont;
// 
// 	max = gsl_vector_max(eval);
// 	min = gsl_vector_min(eval);
// 	
// 	interval = fabs(max - min);
// 	points = interval/(1.0*M);
// 	
// 	for ( i = 0 ; i < M ; i++ )
// 	{
// 		val1 = min + i*points;
// 		val2 = min + (i+1)*points;
// 		valm = (val1 + val2)/2.0;
// 		cont = 0;
// 		for ( j = 0 ; j < n ; j++ )
// 		{
// 			if ( (gsl_vector_get(eval, j) >= val1) && (gsl_vector_get(eval, j) < val2))
// 			{
// 				cont++;
// 			}
// 		}
// 		gsl_matrix_set(hist, i, 0, valm);
// 		gsl_matrix_set(hist, i, 1, cont);
// 	}	
// }

/** gcc -lm -lgsl -O3 -lgslcblas autovalores-matriz.c -o autovalores-matriz.exe */